Introduction to Photometric Stereo¶

Read in the Beethoven Data, you will be working with primarily 3 data artifacts:

3D Array: Images (m, n, k)¶

2D binary mask : (m,n)¶

array S : light vectors¶

In [2]:

#read in the file which has I,M,S where I is the images, M is the mask, and S is the light vectors.
V=ps.read_data_file('Beethoven')

Step One: Get aquainted with the data¶

Images, mask, light vectors¶

In [7]:

# Take the I,M,S and store them into variables
Iv=V[0]
mv= V[1]
sv= V[2]

plt.figure(figsize=(7,7))

plt.imshow(Iv)
plt.show()

plt.figure(figsize=(7,7))
plt.imshow(mv)
plt.show()

plt.figure(figsize=(7,7))
plt.imshow(sv)
plt.show()

In [13]:

#displaying the images
for i in range (0,3):
    plt.figure(figsize=(7,7))
    plt.imshow(Iv[:,:,i])
    plt.show()

Step 2: Applying Mask to Images¶

In [16]:

# Take each individual image and apply the mask to it, then reshape the outupt into a (3,nz) array. 
a=Iv[:,:,0]
b=Iv[:,:,1]
c=Iv[:,:,2]

ai=a[nz]
bi=b[nz]
ci=c[nz]

J=np.row_stack((ai,bi,ci))
J.shape

Out[16]:

(3, 28898)

Step 3: Applying Lamberts Equation¶

I is the intensity of each pixel in the image, p is the albedo of each pixel, N is the surface normal and S is the lighting source with the assumption that it is constant over each image.

Take the dot product of the Image times the inverse of the light vector to obtain the Modulated Normal Field.

Step 4: Calculate the Albedo Modulated Normalized Field¶

Step 5: Calculate the Albedo¶

Albedo: material light absorption property

In [39]:

#Store the albedo back into an empty numpy array and visualize the output. 
Beto=np.zeros((256,256))
Beto[nz]=albedo
fig, ax = plt.subplots(figsize=(7,7))
plt.imshow(Beto)
plt.title("Albedo for Beethoven in 2D")
plt.axis('off')

Out[39]:

(-0.5, 255.5, 255.5, -0.5)

In [40]:

#Albedo with a different color scheme
plt.figure(figsize=(7,7))
plt.imshow(Beto, cmap='binary')

Out[40]:

<matplotlib.image.AxesImage at 0x10ef03a90>

Calculate the normal fields¶

Display Normal Field¶

In [46]:

#display images 
plt.figure(figsize=(6,6))
plt.imshow(be1)
plt.title("Normal Field 1")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be2)
plt.title("Normal Field 2")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be3)
plt.title("Normal Field 3")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(Betov)
plt.title("All Normals")
plt.show()

In [47]:

#display images 
plt.figure(figsize=(6,6))
plt.imshow(be1, cmap='binary')
plt.title("Normal Field 1")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be2, cmap= 'binary')
plt.title("Normal Field 2")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(be3, cmap='binary')
plt.title("Normal Field 3")
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(Betov, cmap='binary')
plt.title("All Normals")
plt.show()

Step 8: Retrieve 3D¶

In [48]:

#calculate the z to plot using unbiased integrate 
z0=ps.unbiased_integrate(be1,be2,be3,mv)
z0.shape

Out[48]:

(256, 256)

In [49]:

ps.display_depth_matplotlib(z0)

/Users/ViVeri/anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning:

elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison

In [88]:

import plotly.plotly as py
from plotly.graph_objs import *

In [171]:

#the graph is interactive and you can zoom in and out of the 3D figure
plotly.offline.iplot(fig)